###########################################################################################
# Replication Code for 
# Helbling, Maxwell & Traunmüller: 
# Numbers, Selectivity and Rights. The Conditional Nature of Immigration Policy Preferences
# Comparative Political Studies
# 24.03.2023
# Contact: traunmueller@uni-mannheim.de 
###########################################################################################

library(arm)
library(readstata13)
library(dplyr)
library(multiwayvcov)
library(lmtest)
library(RColorBrewer)
library(xtable)
library(stargazer)
library(blme)

Sys.setlocale("LC_CTYPE", "UTF-8")
setwd("~/YOURPATH/")
data <- read.dta13("replication_data_numbers_selectivity_rights_study_1_2.dta")

###########################################################################################
# Recode Variables
###########################################################################################

        data$agree_vig[data$agree_vig==0] <- NA
        data$agree_vig[data$agree_vig==9] <- NA
        data$agree_vig <- data$agree_vig*(-1) + 6

        data$agree_vig.2 <- ifelse(data$agree_vig>=4, 1, 0)

        #First policy item: high immigration number=1   
        data$highnumb <- NA
        data$highnumb <- ifelse(data$Gruppe_1==1 | data$Gruppe_1==2 | data$Gruppe_1==3 | data$Gruppe_1==4, 1, 0)

        #Second policy item: liberal immigration requirements=1
        data$easyacc <- NA
        data$easyacc <- ifelse(data$Gruppe_1==1 | data$Gruppe_1==3 | data$Gruppe_1==5 | data$Gruppe_1==8, 1, 0)

        #Third policy item: generous migrant rights=1
        data$morerights <- NA
        data$morerights <- ifelse(data$Gruppe_1==1 | data$Gruppe_1==2 | data$Gruppe_1==5 | data$Gruppe_1==7, 1, 0)

        # Full Policy Interaction
        data$int.pol <- interaction(data$highnumb, data$easyacc, data$morerights)
        levels(data$int.pol) <- c("RRR", "LRR", "RLR", "LLR", "RRL", "LRL", "RLL", "LLL")

        # Pro-Immigrant Subset
        data$pro_immig <- data$v_240
        data$pro_immig[data$v_240==0] <- NA
        data$pro_immig <- ifelse(data$pro_immig<=3, 1, data$pro_immig)
        data$pro_immig <- ifelse(data$pro_immig>=6, -1, data$pro_immig)
        data$pro_immig <- ifelse(data$pro_immig==4 | data$pro_immig==5, 0, data$pro_immig)


        # Covariates
        data$r.female <- ifelse(data$v_2==2, 1, 0)
        data$r.age <- (2020-data$v_3)
        data$r.abi  <- ifelse(data$v_4==4, 1, 0)  
        data$r.german <- ifelse(data$v_9==1, 1, 0)
        data$v_311 <- ifelse(data$v_311==-77, NA, data$v_311)
        data$r.nation <- data$v_311
        data$r.immig  <- data$v_240
        data$r.leftright  <- data$v_121

############################################################################################
# Descriptive Statistics
############################################################################################
        
#**********
# Table A1:
        stargazer(data[data$wave<=13,c("r.female", "r.age", "r.abi", "r.german", "r.nation", "r.immig", "r.leftright")])

#**********
# Table A3:
        stargazer(data[data$wave>=14,c("r.female", "r.age", "r.abi", "r.german", "r.nation", "r.immig", "r.leftright")])

############################################################################################
# Balance Testing
############################################################################################
        
        b.1 <- lm(r.female ~ int.pol, data=data[data$wave<=13,])
          summary(b.1)
          anova(b.1)
        b.2 <- lm(r.age ~ int.pol, data=data[data$wave<=13,])
          summary(b.2)
          anova(b.2)
        b.3 <- lm(r.abi ~ int.pol, data=data[data$wave<=13,])
          summary(b.3)
          anova(b.3)
        b.4 <- lm(r.german ~ int.pol, data=data[data$wave<=13,])
          summary(b.4)
          anova(b.4)
        b.5 <- lm(r.leftright ~ int.pol, data=data[data$wave<=13,])
          summary(b.5)
          anova(b.5)
        b.6 <- lm(r.nation ~ int.pol, data=data[data$wave<=13,])
          summary(b.6)
          anova(b.6)
        b.7 <- lm(r.immig ~ int.pol, data=data[data$wave<=13,])
          summary(b.7)
          anova(b.7)

        rbind(anova(b.1), anova(b.2), anova(b.3), anova(b.4), 
        anova(b.5), anova(b.6), anova(b.7))

        balance.test <- rbind(anova(b.1)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.2)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.3)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.4)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.5)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.6)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.7)[1,c("Df", "F value", "Pr(>F)")])
        rownames(balance.test) <- c("Female", "Age", "Education", "German", "Left-Right-Ideology", "National Attachment", "Immigration Attitude")

#**********
# Table A5:
        xtable(balance.test)


        b.1 <- lm(r.female ~ int.pol, data=data[data$wave>=14,])
          summary(b.1)
          anova(b.1)
        b.2 <- lm(r.age ~ int.pol, data=data[data$wave>=14,])
          summary(b.2)
          anova(b.2)
        b.3 <- lm(r.abi ~ int.pol, data=data[data$wave>=14,])
          summary(b.3)
          anova(b.3)
        b.4 <- lm(r.german ~ int.pol, data=data[data$wave>=14,])
          summary(b.4)
          anova(b.4)
        b.5 <- lm(r.leftright ~ int.pol, data=data[data$wave>=14,])
          summary(b.5)
          anova(b.5)
        b.6 <- lm(r.nation ~ int.pol, data=data[data$wave>=14,])
          summary(b.6)
          anova(b.6)
        b.7 <- lm(r.immig ~ int.pol, data=data[data$wave>=14,])
          summary(b.7)
          anova(b.7)

        rbind(anova(b.1), anova(b.2), anova(b.3), anova(b.4), 
        anova(b.5), anova(b.6), anova(b.7))

        balance.test <- rbind(anova(b.1)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.2)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.3)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.4)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.5)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.6)[1,c("Df", "F value", "Pr(>F)")],
                      anova(b.7)[1,c("Df", "F value", "Pr(>F)")])
        rownames(balance.test) <- c("Female", "Age", "Education", "German", "Left-Right-Ideology", "National Attachment", "Immigration Attitude")

#**********
# Table A6:
        xtable(balance.test)


############################################################################################
# Prefences
###########################################################################################
        
        table(data$v_442)   # Borders

        table(data$v_443)   # Equal Rights


############################################################################################
# Simple Description
############################################################################################
        
      raw.means <- aggregate(agree_vig ~ int.pol, data=data[data$wave<=13,], mean)

      mean(data$agree_vig, na.rm=T)
      sd(data$agree_vig, na.rm=T)

      raw.means.2 <- aggregate(agree_vig.2 ~ int.pol, data=data[data$wave<=13,], mean)

      raw.means.l <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==1 & data$wave<=13,], mean)
      raw.means.c <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==0 & data$wave<=13,], mean)
      raw.means.r <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==-1 & data$wave<=13,], mean)

      agree.dat <- cbind(raw.means.l, raw.means.c[,2], raw.means.r[,2])
      colnames(agree.dat) <- c("policy", "pro", "neutral", "anti")
      agree.dat <- as.matrix(agree.dat[,2:4])

      se.raw.means.l <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==1 & data$wave<=13,], sd)
      se.raw.means.c <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==0 & data$wave<=13,], sd)
      se.raw.means.r <- aggregate(agree_vig.2 ~ int.pol, data=data[data$pro_immig==-1 & data$wave<=13,], sd)

      n.raw.means.l <- table(data$int.pol[data$pro_immig==1 & data$wave<=13 & is.na(data$agree_vig.2)==F])
      n.raw.means.c <- table(data$int.pol[data$pro_immig==0 & data$wave<=13 & is.na(data$agree_vig.2)==F])
      n.raw.means.r <- table(data$int.pol[data$pro_immig==-1 & data$wave<=13 & is.na(data$agree_vig.2)==F])

      se.agree.dat <- cbind(se.raw.means.l[,1], se.raw.means.l[,2]/sqrt(n.raw.means.l), se.raw.means.c[,2]/sqrt(n.raw.means.c), se.raw.means.r[,2]/sqrt(n.raw.means.r))
      colnames(se.agree.dat) <- c("policy", "pro", "neutral", "anti")
      se.agree.dat <- as.matrix(se.agree.dat[,2:4])

######################################################################################################
# STUDY 1
######################################################################################################
# Plot Relative Distributions Across Treatment Groups.
######################################################################################################
      
      policy.levels <- c("Tighten requirements \n Restrict rights", "Tighten requirements \n Restrict rights", 
                   "Lower requirements \n Restrict rights", "Lower requirements \n Restrict rights", 
                   "Tighten requirements \n Extend rights", "Tighten requirements \n Extend rights", 
                   "Lower requirements \n Extend rights", "Lower requirements \n Extend rights")
#**********
# Figure 1
      par(mgp=c(2, .2, 0), mfrow=c(2,4), mar=c(4,2,4,2), oma=c(2,2,3,2), cex.main=.8, cex.axis=.8, las=1, yaxs="r")

      for(i in c(8, 4, 6, 2, 7, 3, 5, 1)){
          p <- barplot(agree.dat[i, 1:3], main=policy.levels[i], ylim=c(0, .75), axes=F, col=rev(brewer.pal(3, "PiYG")), border="darkgrey")
          axis(2, col="white")
          #text(p, -0.03, round(agree.dat[i, 1:3], 2), pos=3, cex=.8, col=c("black", "darkgrey", "black"))
          grid(lty=1, lwd=.5, col="lightgrey")
          arrows(p, agree.dat[i, 1:3]-1.96*se.agree.dat[i, 1:3], p, agree.dat[i, 1:3]+1.96*se.agree.dat[i, 1:3], code=3, angle=90, length=.05, col="darkgrey")
        }

      mtext("A) Increase Numbers", side=3, line=-0.5, adj=0, outer=T, font=1, cex=.8)
      mtext("B) Decrease Numbers", side=3, line=-21, adj=0, outer=T, font=1, cex=.8)


#########################################################################################
# Regression Models
######################################################################################################
      
      m.0 <- lm(agree_vig ~ highnumb + easyacc + morerights + factor(wave), data=data[data$wave<=13,]) 
      summary(m.0)

      vcov_wave <- cluster.vcov(m.0, cbind(data$wave[data$wave<=13]))
      coeftest(m.0 , vcov_wave)

      m.1 <- lm(agree_vig ~ highnumb + easyacc + morerights + highnumb:easyacc + highnumb:morerights +  factor(wave), data=data[data$wave<=13,]) 
      summary(m.1)

      m.2 <- lm(agree_vig ~ highnumb + easyacc + morerights + highnumb:easyacc + highnumb:morerights + easyacc:morerights + factor(wave), data=data[data$wave<=13,]) 
      summary(m.2)
      anova(m.0, m.1, m.2)


      # Print Regression Table
      stargazer(m.0, m.1, m.2, omit.stat=c("ser","adj.rsq", "f"), 
          no.space=TRUE,
          keep = c("Constant", "highnumb", "easyacc", "morerights"),
          intercept.bottom=T,
          covariate.labels = c("Increase number", "Lower requirements", "Extend rights", "Increase number*Lower requirements", "Increase number*Extend rights", "Lower requirements*Extend rights", "Increase number*Lower requirements*Extend rights"),
          dep.var.labels=c("Agreement with Immigration Policy \n (Scale: 1-5, Mean: 2.73, SD: 1.21)"),
          digits=2)

      # Simulation
      sim.m <- sim(m.1, 5000)

      # Combine terms
      LLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + 
        coef(sim.m)[,"highnumb:easyacc"] + coef(sim.m)[,"highnumb:morerights"] + coef(sim.m)[,"easyacc:morerights"] + 
        coef(sim.m)[,"highnumb:easyacc:morerights"]

      LLR <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"highnumb:easyacc"] 

      LRL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"highnumb:morerights"] 

      LRR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"highnumb"] 

      RLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"easyacc:morerights"] 

      RLR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"easyacc"] 

      RRL <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"morerights"] 

      RRR <- coef(sim.m)[,"(Intercept)"]
  
      result.mat <- cbind(LLL, LLR, LRL, LRR, RLL, RLR, RRL, RRR)  
  
      result.means <- apply(result.mat, 2, mean)
      result.ci.up <- apply(result.mat, 2, quantile, .975)
      result.ci.lo <- apply(result.mat, 2, quantile, .025)
      result.ci.up.83 <- apply(result.mat, 2, quantile, .915)
      result.ci.lo.83 <- apply(result.mat, 2, quantile, .085)


      pol.sel <- c("lower requirements", "lower requirements", "tighten requirements", "tighten requirements")
      pol.rig <- c("extend rights", "restrict rights", "extend rights", "restrict rights")

      par(las=1, mgp=c(2, .5, 0), mar=c(4, 3, 3, 3), mfrow=c(1,1))
      plot(result.means[5:8], type="o", pch=19, ylab="Agreement with Immigration Policy (1-5)", cex.lab=.8, xlab="", axes=F, ylim=c(2, 3.5), col=brewer.pal(5, "PiYG")[1])
      axis(2, col="white", cex.axis=.8)
      axis(4, col="white", cex.axis=.8)
      grid(lty=1, lwd=.5, col="lightgrey")
      axis(1, at=1:4, labels=pol.sel, col="white", cex.axis=.8, line=.5)
      axis(1, at=1:4, labels=pol.rig, col="white", cex.axis=.8, line=1)
      points(result.means[1:4], type="o", pch=19, col=brewer.pal(5, "PiYG")[5])
      text(2.5, 3.1, "decrease numbers", col=brewer.pal(5, "PiYG")[1], cex=.8)
      text(2.5, 2.25, "increase numbers", col=brewer.pal(5, "PiYG")[5], cex=.8)

      segments(1:4, result.ci.lo[5:8], 1:4, result.ci.up[5:8], col=brewer.pal(5, "PiYG")[1])
      segments(1:4, result.ci.lo[1:4], 1:4, result.ci.up[1:4], col=brewer.pal(5, "PiYG")[5])
      segments(1:4, result.ci.lo.83[5:8], 1:4, result.ci.up.83[5:8], col=brewer.pal(5, "PiYG")[1], lwd=2)
      segments(1:4, result.ci.lo.83[1:4], 1:4, result.ci.up.83[1:4], col=brewer.pal(5, "PiYG")[5], lwd=2)

###########################################################################################
# Multiple Comparison Plot
###########################################################################################

#**********      
# Figure A8:
      policy.levels.2 <- c("lower requirements \n extend rights",
                     "lower requirements \n restrict rights",
                     "tighten requirements \n extend rights",
                     "tighten requirements \n restrict rights",
                     "lower requirements \n extend rights",
                     "lower requirements \n restrict rights",
                     "tighten requirements \n extend rights",
                     "tighten requirements \n restrict rights"
                     )

      par(mgp=c(2, 1, 0), fg="black", mfrow=c(1,1), cex.axis=.6,  oma=c(3,6,0,0))
      multicomp.plot(result.mat, alpha=.05, label=c("", "", "", "", "", "", "", ""), show.pvalue=T, vertical.line=F, label.on.which.axis = 3, col.same="#f7f7f7", col.high="#b8e186", col.low="#f1b6da", label.as.shortlabel=TRUE, main="", mar=c(2,3,2,2))
      axis(2, at=1:8, labels=rev(policy.levels.2), tck=0, line=-.5, col=NA)
      axis(1, at=1:8, labels=(policy.levels.2), tck=0, line=-.5, col=NA, las=2)

      rect(0.5, 4.5, 4.5, 10.5, xpd=NA)
      rect(4.5, 4.5, 8.5, 10.5, xpd=NA)
      text(2.5, 10.25, "Increase Numbers", xpd=NA, pos=1, cex=.8)
      text(6.5, 10.25, "Decrease Numbers", xpd=NA, pos=1, cex=.8)

      rect(-1.5, 4.5, 4.5, 8.5, xpd=NA)
      rect(-1.5, 0.5, 4.5, 4.5, xpd=NA)
      text(-1, 7, "Increase \nNumbers", xpd=NA, pos=1, cex=.8, padj=2)
      text(-1, 3, "Decrease\nNumbers", xpd=NA, pos=1, cex=.8)

      mtext("A) Pro-Immigration", side=3, line=4.5, at=-.5)
      mtext("B) Neutral", side=3, line=4.5, at=-.3)
      mtext("C) Anti-Immigration", side=3, line=4.5, at=-.35)

############################################################################################
# Bayesian Hierarchical Models
############################################################################################
      
      mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$wave<=13,])
      summary(mlm.0)

      mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$wave<=13,])
      summary(mlm.1)
      anova(mlm.0, mlm.1)

      # Pro-Immigration
      mlm.0.l <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$pro_immig==1 & data$wave<=13,])
      summary(mlm.0.l)

      mlm.1.l <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$pro_immig==1 & data$wave<=13,])
      summary(mlm.1.l)
      anova(mlm.0.l, mlm.1.l)

      # Neutral
      mlm.0.n <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$pro_immig==0 & data$wave<=13,])
      summary(mlm.0.n)

      mlm.1.n <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$pro_immig==0 & data$wave<=13,])
      summary(mlm.1.n)
      anova(mlm.0.n, mlm.1.n)


      # Anti-Immigration
      mlm.0.r <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$pro_immig==-1 & data$wave<=13,])
      summary(mlm.0.r)

      mlm.1.r <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$pro_immig==-1 & data$wave<=13,])
      summary(mlm.1.r)
      anova(mlm.0.r, mlm.1.r)

#*********
# Table 1
      # Model Comparison Table
      xtable(rbind(anova(mlm.0, mlm.1)[,5:8],  
      anova(mlm.0.l, mlm.1.l)[,5:8], 
      anova(mlm.0.r, mlm.1.r)[,5:8]))

#*********     
# Figure 2
      par(mfrow=c(1,4), mar=c(3, 1, 1, 1))
      names <- c("Increase Numbers", "Lower Requirements", "Extend Rights", "", "SD: Full Policy Interactions", "SD: Waves", "SD: Residuals")

      plot(rep(1,7), 7:1, type="n", ann=F, axes=F)

      # Full Sample
      est.0 <- c(fixef(mlm.0)[-1], NA, NA,
           sqrt(summary(mlm.0)$varcor$wave[1]),
           summary(mlm.0)$sigma)

      est.1 <- c(fixef(mlm.1)[-1], NA,
           sqrt(summary(mlm.1)$varcor$int.pol[1]),
           sqrt(summary(mlm.1)$varcor$wave[1]),
           summary(mlm.1)$sigma)

      plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Full Sample")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0)[-1]-1.96*se.fixef(mlm.0)[-1], c(7:5)-.2, fixef(mlm.0)[-1]+1.96*se.fixef(mlm.0)[-1], c(7:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      axis(2, at=7:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
      segments(fixef(mlm.1)[-1]-1.96*se.fixef(mlm.1)[-1], c(7:5)+.2, fixef(mlm.1)[-1]+1.96*se.fixef(mlm.1)[-1], c(7:5)+.2)
      text(est.1[5], 3+.2, "***",  pos=4)

      # Pro Immigration
      est.0 <- c(fixef(mlm.0.l)[-1], NA, NA,
           sqrt(summary(mlm.0.l)$varcor$wave[1]),
           summary(mlm.0.l)$sigma)

      est.1 <- c(fixef(mlm.1.l)[-1], NA,
           sqrt(summary(mlm.1.l)$varcor$int.pol[1]),
           sqrt(summary(mlm.1.l)$varcor$wave[1]),
           summary(mlm.1.l)$sigma)

      plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Pro-Immigration")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0.l)[-1]-1.96*se.fixef(mlm.0.l)[-1], c(7:5)-.2, fixef(mlm.0.l)[-1]+1.96*se.fixef(mlm.0.l)[-1], c(7:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      #axis(2, at=9:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
      segments(fixef(mlm.1.l)[-1]-1.96*se.fixef(mlm.1.l)[-1], c(7:5)+.2, fixef(mlm.1.l)[-1]+1.96*se.fixef(mlm.1.l)[-1], c(7:5)+.2)


      # Anti-Immigration
      est.0 <- c(fixef(mlm.0.r)[-1], NA, NA,
           sqrt(summary(mlm.0.r)$varcor$wave[1]),
           summary(mlm.0.r)$sigma)

      est.1 <- c(fixef(mlm.1.r)[-1], NA,
           sqrt(summary(mlm.1.r)$varcor$int.pol[1]),
           sqrt(summary(mlm.1.r)$varcor$wave[1]),
           summary(mlm.1.r)$sigma)

      plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Anti-Immigration")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0.r)[-1]-1.96*se.fixef(mlm.0.r)[-1], c(7:5)-.2, fixef(mlm.0.r)[-1]+1.96*se.fixef(mlm.0.r)[-1], c(7:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      #axis(2, at=9:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
      segments(fixef(mlm.1.r)[-1]-1.96*se.fixef(mlm.1.r)[-1], c(7:5)+.2, fixef(mlm.1.r)[-1]+1.96*se.fixef(mlm.1.r)[-1], c(7:5)+.2)
      text(est.1[5], 3+.2, "***",  pos=4)

      # Specific Policies:

      raw.means.2 <- aggregate(agree_vig.2 ~ int.pol, data=data[data$wave<=13,], mean)

      mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$wave<=13,])
      summary(mlm.1)

      coef(mlm.1)$int.pol

      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RRR",]), na.rm=T) + coef(mlm.1)$int.pol["RRR",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LRR",]), na.rm=T) + coef(mlm.1)$int.pol["LRR",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RLR",]), na.rm=T) + coef(mlm.1)$int.pol["RLR",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LLR",]), na.rm=T) + coef(mlm.1)$int.pol["LLR",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RRL",]), na.rm=T) + coef(mlm.1)$int.pol["RRL",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LRL",]), na.rm=T) + coef(mlm.1)$int.pol["LRL",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RLL",]), na.rm=T) + coef(mlm.1)$int.pol["RLL",1]
      mean(fixef(mlm.1)[-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LLL",]), na.rm=T) + coef(mlm.1)$int.pol["LLL",1]

      s.mlm.1 <- sim(mlm.1, 1000)
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RRR",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"RRR",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LRR",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"LRR",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RLR",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"RLR",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LLR",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"LLR",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RRL",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"RRL",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LRL",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"LRL",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="RLL",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"RLL",] + fixef(s.mlm.1)[,1]
      apply(fixef(s.mlm.1)[,-1]%*%t(cbind(data$highnumb, data$easyacc, data$morerights)[data$int.pol=="LLL",]), 1, mean, na.rm=T) + coef(s.mlm.1)$ranef$int.pol[,"LLL",] + fixef(s.mlm.1)[,1]

######################################################################################
# ROBUSTNESS
######################################################################################
      
    # Left
    mlm.0.l <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$r.leftright<6 & data$wave<=13,])
    summary(mlm.0.l)

    mlm.1.l <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$r.leftright<6 & data$wave<=13,])
    summary(mlm.1.l)
    anova(mlm.0.l, mlm.1.l)

    # Middle
    mlm.0.n <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$r.leftright==6 & data$wave<=13,])
    summary(mlm.0.n)

    mlm.1.n <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$r.leftright==6 & data$wave<=13,])
    summary(mlm.1.n)
    anova(mlm.0.n, mlm.1.n)


    #  Right
    mlm.0.r <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$r.leftright>6 & data$wave<=13,])
    summary(mlm.0.r)

    mlm.1.r <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$r.leftright>6  & data$wave<=13,])
    summary(mlm.1.r)
    anova(mlm.0.r, mlm.1.r)

#**********
# Table A9:
    # Model Comparison Table
    xtable(rbind(
             anova(mlm.0.l, mlm.1.l)[,5:8], 
             anova(mlm.0.n, mlm.1.n)[,5:8],  
             anova(mlm.0.r, mlm.1.r)[,5:8]))


#***********
# Figure A6:
    par(mfrow=c(1,4), mar=c(3, 1, 1, 1))
    names <- c("Increase Numbers", "Lower Requirements", "Extend Rights", "", "SD: Full Policy Interactions", "SD: Waves", "SD: Residuals")

    plot(rep(1,7), 7:1, type="n", ann=F, axes=F)

    # Left
    est.0 <- c(fixef(mlm.0.l)[-1], NA, NA,
           sqrt(summary(mlm.0.l)$varcor$wave[1]),
           summary(mlm.0.l)$sigma)

    est.1 <- c(fixef(mlm.1.l)[-1], NA,
           sqrt(summary(mlm.1.l)$varcor$int.pol[1]),
           sqrt(summary(mlm.1.l)$varcor$wave[1]),
           summary(mlm.1.l)$sigma)

    plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Left")
    grid(lty=1, lwd=.5, col="lightgrey")
    segments(fixef(mlm.0.l)[-1]-1.96*se.fixef(mlm.0.l)[-1], c(7:5)-.2, fixef(mlm.0.l)[-1]+1.96*se.fixef(mlm.0.l)[-1], c(7:5)-.2)
    abline(v=0, lty=1, lwd=1.5, col="grey")
    axis(2, at=7:1, names, lty=0)
    axis(1, lty=0)
    points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
    points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
    segments(fixef(mlm.1.l)[-1]-1.96*se.fixef(mlm.1.l)[-1], c(7:5)+.2, fixef(mlm.1.l)[-1]+1.96*se.fixef(mlm.1.l)[-1], c(7:5)+.2)
    text(est.1[5], 3+.2, "***",  pos=4)

    # Neutral
    est.0 <- c(fixef(mlm.0.n)[-1], NA, NA,
           sqrt(summary(mlm.0.n)$varcor$wave[1]),
           summary(mlm.0.n)$sigma)

    est.1 <- c(fixef(mlm.1.n)[-1], NA,
           sqrt(summary(mlm.1.n)$varcor$int.pol[1]),
           sqrt(summary(mlm.1.n)$varcor$wave[1]),
           summary(mlm.1.n)$sigma)

    plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Middle")
    grid(lty=1, lwd=.5, col="lightgrey")
    segments(fixef(mlm.0.n)[-1]-1.96*se.fixef(mlm.0.n)[-1], c(7:5)-.2, fixef(mlm.0.n)[-1]+1.96*se.fixef(mlm.0.n)[-1], c(7:5)-.2)
    abline(v=0, lty=1, lwd=1.5, col="grey")
    #axis(2, at=7:1, names, lty=0)
    axis(1, lty=0)
    points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
    points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
    segments(fixef(mlm.1.n)[-1]-1.96*se.fixef(mlm.1.n)[-1], c(7:5)+.2, fixef(mlm.1.n)[-1]+1.96*se.fixef(mlm.1.n)[-1], c(7:5)+.2)
    text(est.1[5], 3+.2, "***",  pos=4)

    # Right
    est.0 <- c(fixef(mlm.0.r)[-1], NA, NA,
           sqrt(summary(mlm.0.r)$varcor$wave[1]),
           summary(mlm.0.r)$sigma)

    est.1 <- c(fixef(mlm.1.r)[-1], NA,
           sqrt(summary(mlm.1.r)$varcor$int.pol[1]),
           sqrt(summary(mlm.1.r)$varcor$wave[1]),
           summary(mlm.1.r)$sigma)

    plot(est.0, c(7:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 7.5), xlim=c(-.5, .5), main="Right")
    grid(lty=1, lwd=.5, col="lightgrey")
    segments(fixef(mlm.0.r)[-1]-1.96*se.fixef(mlm.0.r)[-1], c(7:5)-.2, fixef(mlm.0.r)[-1]+1.96*se.fixef(mlm.0.r)[-1], c(7:5)-.2)
    abline(v=0, lty=1, lwd=1.5, col="grey")
    #axis(2, at=9:1, names, lty=0)
    axis(1, lty=0)
    points(est.0, c(7:1)-.2, pch=19, cex=.6, col="white")
    points(est.1, c(7:1)+.2, pch=19, col=c(rep(1, 4), "maroon3", rep(1,2)))
    segments(fixef(mlm.1.r)[-1]-1.96*se.fixef(mlm.1.r)[-1], c(7:5)+.2, fixef(mlm.1.r)[-1]+1.96*se.fixef(mlm.1.r)[-1], c(7:5)+.2)
    text(est.1[5], 3+.2, "***",  pos=4)

######################################################################################
# ROBUSTNESS:
# Different Cut-Offs for Anti and Pro Immigration
######################################################################################
    
    # Pro-Immigration
    est.mat.l <- matrix(NA, 4, 6)
    se.mat.l <- matrix(NA, 4, 3)
    p.mat.l <- matrix(NA, 4, 1)
    coef.mat.l <- matrix(NA, 4, 8)

    for(i in 1:4){
  
      data$pro_immig2 <- data$v_240
      data$pro_immig2[data$v_240==0] <- NA
      data$pro_immig2 <- ifelse(data$pro_immig2<=i, 1, data$pro_immig2)
      data$pro_immig2 <- ifelse(data$pro_immig2>=9-i, -1, data$pro_immig2)
  
  
      mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$pro_immig2==1 & data$wave<=13,])
      summary(mlm.0)
  
      mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$pro_immig2==1 & data$wave<=13,])
      summary(mlm.1.l)
      anova(mlm.0, mlm.1)
  
  
      est.mat.l[i,]  <- c(fixef(mlm.1)[-1],
                    sqrt(summary(mlm.1)$varcor$int.pol[1]),
                    sqrt(summary(mlm.1)$varcor$wave[1]),
                    summary(mlm.1)$sigma)
      se.mat.l[i,]  <- se.fixef(mlm.1)[-1]
  
      p.mat.l[i,] <- anova(mlm.0, mlm.1)$Pr[2]
      coef.mat.l[i,] <- coef(mlm.1)$int.pol[,1]
      }


      # Anti-Immigration
      est.mat.r <- matrix(NA, 4, 6)
      se.mat.r <- matrix(NA, 4, 3)
      p.mat.r <- matrix(NA, 4, 1)

      for(i in 1:4){
  
        data$pro_immig2 <- data$v_240
        data$pro_immig2[data$v_240==0] <- NA
        data$pro_immig2 <- ifelse(data$pro_immig2<=i, 1, data$pro_immig2)
        data$pro_immig2 <- ifelse(data$pro_immig2>=9-i, -1, data$pro_immig2)
  
  
      mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$pro_immig2==-1 & data$wave<=13,])
      summary(mlm.0)
  
      mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$pro_immig2==-1 & data$wave<=13,])
      summary(mlm.1)
      anova(mlm.0, mlm.1)
  
  
      est.mat.r[i,]  <- c(fixef(mlm.1)[-1],
                    sqrt(summary(mlm.1)$varcor$int.pol[1]),
                    sqrt(summary(mlm.1)$varcor$wave[1]),
                    summary(mlm.1)$sigma)
      se.mat.r[i,]  <- se.fixef(mlm.1)[-1]
  
      p.mat.r[i,] <- anova(mlm.0, mlm.1)$Pr[2]
  
      }

#***********
# Figure A4:

      par(mfrow=c(2,3), mar=c(3, 3, 1, 1), family="Gill Sans", font.main=1, cex.main=1, cex.axis=.8, cex.lab=.8)

      plot(1:9, c(est.mat.l[,1], NA, rev(est.mat.r[,1])), pch=19, ylim=c(-.4, .4), main="Increase Numbers", axes=F, xlab="Scale Cutoff Points", ylab="Est.")
      axis(1, at=1:9, labels=c("<=0", "<=1", "<=2", "<=3", "", ">=4", ">=5", ">=6", ">=7"))
      axis(2)
      segments(1:9, c(est.mat.l[,1]-1.96*se.mat.l[,1], NA, rev(est.mat.r[,1]-1.96*se.mat.r[,1])), 1:9, c(est.mat.l[,1]+1.96*se.mat.l[,1], NA, rev(est.mat.r[,1]+1.96*se.mat.r[,1])))
      abline(h=0, col="grey")
      text(2.5, 0, "Pro-Immigration", cex=.8, pos=1)
      text(7.5, 0, "Anti-Immigration", cex=.8, pos=3)

      plot(1:9, c(est.mat.l[,2], NA, rev(est.mat.r[,2])), pch=19, ylim=c(-.4, .4), main="Lower Requirements", axes=F, xlab="Scale Cutoff Points", ylab="Est.")
      axis(1, at=1:9, labels=c("<=0", "<=1", "<=2", "<=3", "", ">=4", ">=5", ">=6", ">=7"))
      axis(2)
      segments(1:9, c(est.mat.l[,2]-1.96*se.mat.l[,2], NA, rev(est.mat.r[,2]-1.96*se.mat.r[,2])), 1:9, c(est.mat.l[,2]+1.96*se.mat.l[,2], NA, rev(est.mat.r[,2]+1.96*se.mat.r[,2])))
      abline(h=0, col="grey")
      text(2.5, 0, "Pro-Immigration", cex=.8, pos=1)
      text(7.5, 0, "Anti-Immigration", cex=.8, pos=3)

      plot(1:9, c(est.mat.l[,3], NA, rev(est.mat.r[,3])), pch=19, ylim=c(-.4, .4), main="Extend Rights", axes=F, xlab="Scale Cutoff Points", ylab="Est.")
      axis(1, at=1:9, labels=c("<=0", "<=1", "<=2", "<=3", "", ">=4", ">=5", ">=6", ">=7"))
      axis(2)
      segments(1:9, c(est.mat.l[,3]-1.96*se.mat.l[,3], NA, rev(est.mat.r[,3]-1.96*se.mat.r[,3])), 1:9, c(est.mat.l[,3]+1.96*se.mat.l[,3], NA, rev(est.mat.r[,3]+1.96*se.mat.r[,3])))
      abline(h=0, col="grey")
      text(2.5, 0, "Pro-Immigration", cex=.8, pos=1)
      text(7.5, 0, "Anti-Immigration", cex=.8, pos=3)

      plot(1:9, c(est.mat.l[,4], NA, rev(est.mat.r[,4])), pch=19, ylim=c(0, .2), main="SD: Full Policy Interactions", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
      axis(1, at=1:9, labels=c("<=0", "<=1", "<=2", "<=3", "", ">=4", ">=5", ">=6", ">=7"))
      axis(2)
      abline(h=0, col="grey")
      text(2.5, 0.1, "Pro-Immigration", cex=.8, pos=3)
      text(7.5, 0.1, "Anti-Immigration", cex=.8, pos=1)

      plot(1:9, c(p.mat.l[,1], NA, rev(p.mat.r[,1])), pch=19, ylim=c(0, 1), main="Pr(>Chisq)", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
      axis(1, at=1:9, labels=c("<=0", "<=1", "<=2", "<=3", "", ">=4", ">=5", ">=6", ">=7"))
      axis(2)
      abline(h=0.05, col="grey")
      text(2.5, .5, "Pro-Immigration", cex=.8, pos=1)
      text(7.5, 0, "Anti-Immigration", cex=.8, pos=3)

######################################################################################
# ROBUSTNESS: Wave-Specific: Waves 1-13
######################################################################################
      
      # Full Sample
      est.mat <- matrix(NA, 10, 6)
      se.mat <- matrix(NA, 10, 3)
      p.mat <- matrix(NA, 10, 1)

      for(i in 1:10){
        mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$wave>=i & data$wave<=i+3,])
        summary(mlm.0)

        mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$wave>=i & data$wave<=i+3,])
        summary(mlm.1)
        anova(mlm.0, mlm.1)


        est.mat[i,]  <- c(fixef(mlm.1)[-1],
           sqrt(summary(mlm.1)$varcor$int.pol[1]),
           sqrt(summary(mlm.1)$varcor$wave[1]),
           summary(mlm.1)$sigma)
        se.mat[i,]  <- se.fixef(mlm.1)[-1]

        p.mat[i,] <- anova(mlm.0, mlm.1)$Pr[2]

        }

#***********
# Figure A1:

        par(mfrow=c(2,3), mar=c(3, 3, 1, 1), family="Gill Sans", font.main=1, cex.main=1, cex.axis=.8, cex.lab=.8)

        plot(1:10, est.mat[,1], pch=19, ylim=c(-.4, .2), main="Increase Numbers", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,1]-1.96*se.mat[,1], 1:10, est.mat[,1]+1.96*se.mat[,1])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,2], pch=19, ylim=c(-.4, .2), main="Lower Requirements", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,2]-1.96*se.mat[,2], 1:10, est.mat[,2]+1.96*se.mat[,2])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,3], pch=19, ylim=c(-.4, .2), main="Extend Rights", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,3]-1.96*se.mat[,3], 1:10, est.mat[,3]+1.96*se.mat[,3])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,4], pch=19, ylim=c(0, .2), main="SD: Full Policy Interactions", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0, col="grey")

        plot(1:10, p.mat[,1], pch=19, ylim=c(0, .2), main="Pr(>Chisq)", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0.05, col="grey")

        # Pro-Immigration
       est.mat <- matrix(NA, 10, 6)
        se.mat <- matrix(NA, 10, 3)
        p.mat <- matrix(NA, 10, 1)

        for(i in 1:10){
           mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$wave>=i & data$wave<=i+3 & data$pro_immig==1,])
           summary(mlm.0)
  
           mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$wave>=i & data$wave<=i+3 & data$pro_immig==1,])
           summary(mlm.1)
           anova(mlm.0, mlm.1)
  
  
          est.mat[i,]  <- c(fixef(mlm.1)[-1],
                    sqrt(summary(mlm.1)$varcor$int.pol[1]),
                    sqrt(summary(mlm.1)$varcor$wave[1]),
                    summary(mlm.1)$sigma)
          se.mat[i,]  <- se.fixef(mlm.1)[-1]
  
          p.mat[i,] <- anova(mlm.0, mlm.1)$Pr[2]
       }

#***********
# Figure A2:

        par(mfrow=c(2,3), mar=c(3, 3, 1, 1), family="Gill Sans", font.main=1, cex.main=1, cex.axis=.8, cex.lab=.8)

        plot(1:10, est.mat[,1], pch=19, ylim=c(-.2, .4), main="Increase Numbers", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,1]-1.96*se.mat[,1], 1:10, est.mat[,1]+1.96*se.mat[,1])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,2], pch=19, ylim=c(-.2, .4), main="Lower Requirements", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,2]-1.96*se.mat[,2], 1:10, est.mat[,2]+1.96*se.mat[,2])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,3], pch=19, ylim=c(-.2, .4), main="Extend Rights", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,3]-1.96*se.mat[,3], 1:10, est.mat[,3]+1.96*se.mat[,3])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,4], pch=19, ylim=c(0, .2), main="SD: Full Policy Interactions", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0, col="grey")

        plot(1:10, p.mat[,1], pch=19, ylim=c(0, 1), main="Pr(>Chisq)", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0.05, col="grey")


        # Anti-Immigration
        est.mat <- matrix(NA, 10, 6)
        se.mat <- matrix(NA, 10, 3)
        p.mat <- matrix(NA, 10, 1)

        for(i in 1:10){
             mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$wave>=i & data$wave<=i+3 & data$pro_immig==-1,])
             summary(mlm.0)
  
             mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$wave>=i & data$wave<=i+3 & data$pro_immig==-1,])
              summary(mlm.1)
              anova(mlm.0, mlm.1)
  
  
          est.mat[i,]  <- c(fixef(mlm.1)[-1],
                    sqrt(summary(mlm.1)$varcor$int.pol[1]),
                    sqrt(summary(mlm.1)$varcor$wave[1]),
                    summary(mlm.1)$sigma)
          se.mat[i,]  <- se.fixef(mlm.1)[-1]
  
          p.mat[i,] <- anova(mlm.0, mlm.1)$Pr[2]
  
        }

#***********
# Figure A3:

        par(mfrow=c(2,3), mar=c(3, 3, 1, 1), family="Gill Sans", font.main=1, cex.main=1, cex.axis=.8, cex.lab=.8)

        plot(1:10, est.mat[,1], pch=19, ylim=c(-.4, .2), main="Increase Numbers", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,1]-1.96*se.mat[,1], 1:10, est.mat[,1]+1.96*se.mat[,1])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,2], pch=19, ylim=c(-.4, .2), main="Lower Requirements", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,2]-1.96*se.mat[,2], 1:10, est.mat[,2]+1.96*se.mat[,2])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,3], pch=19, ylim=c(-.4, .2), main="Extend Rights", axes=F, xlab="Waves", ylab="Est.")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        segments(1:10, est.mat[,3]-1.96*se.mat[,3], 1:10, est.mat[,3]+1.96*se.mat[,3])
        abline(h=0, col="grey")

        plot(1:10, est.mat[,4], pch=19, ylim=c(0, .2), main="SD: Full Policy Interactions", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0, col="grey")

        plot(1:10, p.mat[,1], pch=19, ylim=c(0, 1), main="Pr(>Chisq)", axes=F, xlab="Waves", ylab="Est.", col="maroon3")
        axis(1, at=1:10, labels=c("1-4", "2-5", "3-6", "4-7", "5-8", "6-9", "7-10", "8-11", "9-12", "10-13"))
        axis(2)
        abline(h=0.05, col="grey")

######################################################################################
# ROBUSTNESS: Political Parties
######################################################################################
        
#1 - CDU/CSU 
#2 - SPD 
#3 - AfD 
#4 - FDP 
#5 - Die Linke 
#6 - Bündnis 90/Die Grünen 
# Sonstige Parteien 
# Weiß nicht 
# Ich würde nicht wählen gehen 

      est.mat <- matrix(NA, 9, 6)
      se.mat <- matrix(NA, 9, 3)
      p.mat <- matrix(NA, 9, 1)

      for(i in 1:9){
         mlm.0 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | wave), data=data[data$v_122==i & data$wave<=13,])
         summary(mlm.0)
  
         mlm.1 <- blmer(agree_vig.2 ~ highnumb + easyacc + morerights + (1 | int.pol) + (1 | wave), data=data[data$v_122==i & data$wave<=13,])
         summary(mlm.1)
         anova(mlm.0, mlm.1)
  
  
        est.mat[i,]  <- c(fixef(mlm.1)[-1],
                    sqrt(summary(mlm.1)$varcor$int.pol[1]),
                    sqrt(summary(mlm.1)$varcor$wave[1]),
                    summary(mlm.1)$sigma)
        se.mat[i,]  <- se.fixef(mlm.1)[-1]
  
        p.mat[i,] <- anova(mlm.0, mlm.1)$Pr[2]
  
        }

#***********
# Figure A5:

      par(mfrow=c(2,3), mar=c(3, 4, 1, 1), family="Gill Sans", font.main=1, cex.main=1, cex.axis=.8, cex.lab=.8)

      plot( est.mat[,1], 1:9, pch=19, xlim=c(-.4, .2), main="Increase Numbers", axes=F, ylab="", xlab="Est.")
      axis(2, at=1:9, labels=c("CDU/CSU", "SPD", "AfD", "FDP", "Die Linke", "Grüne", "Other", "Don`t know", "Would not vote"))
      axis(1)
      segments(est.mat[,1]-1.96*se.mat[,1], 1:9, est.mat[,1]+1.96*se.mat[,1], 1:9)
      abline(v=0, col="grey")

      plot(est.mat[,2], 1:9,  pch=19, xlim=c(-.4, .2), main="Lower Requirements", axes=F, ylab="", xlab="Est.")
      axis(2, at=1:9, labels=c("CDU/CSU", "SPD", "AfD", "FDP", "Die Linke", "Grüne", "Other", "Don`t know", "Would not vote"))
      axis(1)
      segments(est.mat[,2]-1.96*se.mat[,2], 1:9, est.mat[,2]+1.96*se.mat[,2], 1:9)
      abline(v=0, col="grey")

      plot(est.mat[,3], 1:9,pch=19, xlim=c(-.4, .2), main="Extend Rights", axes=F, ylab="", xlab="Est.")
      axis(2, at=1:9, labels=c("CDU/CSU", "SPD", "AfD", "FDP", "Die Linke", "Grüne", "Other", "Don`t know", "Would not vote"))
      axis(1)
      segments(est.mat[,3]-1.96*se.mat[,3], 1:9, est.mat[,3]+1.96*se.mat[,3], 1:9)
      abline(v=0, col="grey")

      plot(est.mat[,4], 1:9, pch=19, xlim=c(0, .2), main="SD: Full Policy Interactions", axes=F, ylab="", xlab="Est.", col="maroon3")
      axis(2, at=1:9, labels=c("CDU/CSU", "SPD", "AfD", "FDP", "Die Linke", "Grüne", "Other", "Don`t know", "Would not vote"))
      axis(1)
      abline(v=0, col="grey")

      plot(p.mat[,1], 1:9, pch=19, xlim=c(0, 1), main="Pr(>Chisq)", axes=F, ylab="", xlab="Est.", col="maroon3")
      axis(2, at=1:9, labels=c("CDU/CSU", "SPD", "AfD", "FDP", "Die Linke", "Grüne", "Other", "Don`t know", "Would not vote"))
      axis(1)
      abline(v=0.05, col="grey")

###########################################################################################
# Condorcet Voting (Copeland's method)
###########################################################################################
      

      m.2.l <- lm(agree_vig ~ highnumb*easyacc*morerights + factor(wave), data=data[data$pro_immig==1 & data$wave<=13,]) 

      m.2.c <- lm(agree_vig ~ highnumb*easyacc*morerights + factor(wave), data=data[data$pro_immig==0 & data$wave<=13,]) 

      m.2.r <- lm(agree_vig ~ highnumb*easyacc*morerights + factor(wave), data=data[data$pro_immig==-1 & data$wave<=13,]) 


      # Simulation
      sim.m <- sim(m.2.l, 5000)

      # Combine terms
      LLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + 
        coef(sim.m)[,"highnumb:easyacc"] + coef(sim.m)[,"highnumb:morerights"] + coef(sim.m)[,"easyacc:morerights"] + 
        coef(sim.m)[,"highnumb:easyacc:morerights"]

      LLR <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"highnumb:easyacc"] 

      LRL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"highnumb:morerights"] 

      LRR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"highnumb"] 

      RLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"easyacc:morerights"] 

      RLR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"easyacc"] 

      RRL <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"morerights"] 

      RRR <- coef(sim.m)[,"(Intercept)"]

      result.mat <- cbind(LLL, LLR, LRL, LRR, RLL, RLR, RRL, RRR)  


      # Simulation
      sim.m <- sim(m.2.c, 5000)

      # Combine terms
      LLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + 
        coef(sim.m)[,"highnumb:easyacc"] + coef(sim.m)[,"highnumb:morerights"] + coef(sim.m)[,"easyacc:morerights"] + 
        coef(sim.m)[,"highnumb:easyacc:morerights"]

      LLR <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"highnumb:easyacc"] 

      LRL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"highnumb:morerights"] 

      LRR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"highnumb"] 

      RLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"easyacc:morerights"] 

      RLR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"easyacc"] 

      RRL <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"morerights"] 

      RRR <- coef(sim.m)[,"(Intercept)"]

      result.mat <- cbind(LLL, LLR, LRL, LRR, RLL, RLR, RRL, RRR)  

      mat.c <- multicomp.plot(result.mat)

      # Simulation
      sim.m <- sim(m.2.r, 5000)

      # Combine terms
      LLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + 
        coef(sim.m)[,"highnumb:easyacc"] + coef(sim.m)[,"highnumb:morerights"] + coef(sim.m)[,"easyacc:morerights"] + 
        coef(sim.m)[,"highnumb:easyacc:morerights"]

      LLR <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"easyacc"] + coef(sim.m)[,"highnumb:easyacc"] 

      LRL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"highnumb"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"highnumb:morerights"] 

      LRR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"highnumb"] 

      RLL <- coef(sim.m)[,"(Intercept)"] + 
        coef(sim.m)[,"easyacc"] + coef(sim.m)[,"morerights"] + coef(sim.m)[,"easyacc:morerights"] 

      RLR <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"easyacc"] 

      RRL <- coef(sim.m)[,"(Intercept)"] + coef(sim.m)[,"morerights"] 

      RRR <- coef(sim.m)[,"(Intercept)"]

      result.mat <- cbind(LLL, LLR, LRL, LRR, RLL, RLR, RRL, RRR)  

      mat.r <- multicomp.plot(result.mat)


      # Left vs. Right
      condorcet.mat <- mat.l$significant + mat.r$significant 

      condorcet.mat[condorcet.mat>0] <- 1
      condorcet.mat[condorcet.mat<0] <- -1

      condorcet.winner <- apply(condorcet.mat, 2, sum)

#*********
# Table 5:
      xtable(round(t(rbind(condorcet.mat, condorcet.winner))), digits=0)

      condorcet.winner <- sort(condorcet.winner, decreasing=T)
      condorcet.winner



      # Left vs. Right vs. Center
      condorcet.mat <- mat.l$significant + mat.r$significant + mat.c$significant

      condorcet.mat[condorcet.mat>0] <- 1
      condorcet.mat[condorcet.mat<0] <- -1

      condorcet.winner <- apply(condorcet.mat, 2, sum)

      xtable(round(t(rbind(condorcet.mat, condorcet.winner))), digits=0)


      condorcet.winner <- sort(condorcet.winner, decreasing=T)
      condorcet.winner

######################################################################################################
# STUDY 2
######################################################################################################

      data$highnumb2 <- ifelse(data$c_0034=="zu veringern.", 0, 1)

      data$easyacc2 <- ifelse(data$c_0035=="in Bezug auf Bildung versch<e4>rfen" |
                        data$c_0035=="in Bezug auf Bildung verschärfen" |
                        data$c_0035=="in Bezug auf Nationalit<e4>ten versch<e4>rfen" |
                        data$c_0035=="in Bezug auf Nationalitäten verschärfen"
                          , 0, 1)

      # education==0, nationality==1
      data$diffacc <- ifelse(data$c_0035=="in Bezug auf Bildung versch<e4>rfen" |
                          data$c_0035=="in Bezug auf Bildung verschärfen" |
                          data$c_0035=="in Bezug auf Bildung erleichtern" 
                        , 0, 1)

      data$morerights2<- ifelse(data$c_0036=="Rechte einschr<e4>nken, damit sie ihre Sitten und Gebr<e4>uche weniger stark aus<fc>ben k<f6>nnen." |
                          data$c_0036=="Rechte einschr<e4>nken, damit sie weniger stark vom Wohlfahrtsstaat profitieren k<f6>nnen." |
                          data$c_0036=="Rechte einschränken, damit sie ihre Sitten und Gebräuche weniger stark ausüben können." |
                          data$c_0036=="Rechte einschränken, damit sie weniger stark vom Wohlfahrtsstaat profitieren können."
                        , 0, 1)

      # welfare=0, customs=1
      data$diffrights <- ifelse(data$c_0036=="Rechte ausweiten, damit sie st<e4>rker vom Wohlfahrtsstaat profitieren k<f6>nnen." |
                             data$c_0036=="Rechte einschr<e4>nken, damit sie weniger stark vom Wohlfahrtsstaat profitieren k<f6>nnen." |
                             data$c_0036=="Rechte ausweiten, damit sie stärker vom Wohlfahrtsstaat profitieren können." |
                             data$c_0036=="Rechte einschränken, damit sie weniger stark vom Wohlfahrtsstaat profitieren können."
                           , 0, 1)

      data$access[data$easyacc2==1 & data$diffacc==1] <- "lower_nat"
      data$access[data$easyacc2==1 & data$diffacc==0] <- "lower_edu"
      data$access[data$easyacc2==0 & data$diffacc==1] <- "increase_nat"
      data$access[data$easyacc2==0 & data$diffacc==0] <- "increase_edu"

      data$rights[data$morerights2==1 & data$diffrights==1] <- "expand_cust"
      data$rights[data$morerights2==1 & data$diffrights==0] <- "expand_welf"
      data$rights[data$morerights2==0 & data$diffrights==1] <- "restrict_cust"
      data$rights[data$morerights2==0 & data$diffrights==0] <- "restrict_welf"


      data$fullint <- interaction(data$highnumb2, data$access, data$rights)

######################################################################################################
## Models
######################################################################################################
      
      # Full Sample
      mlm.0 <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | wave), data=data[data$wave>13,])
      summary(mlm.0)

      mlm.1<- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | fullint) + (1 | wave), data=data[data$wave>13,])
      summary(mlm.1)
      anova(mlm.0, mlm.1)

      # Pro-Immigration
      mlm.0.l <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | wave), data=data[data$pro_immig==1 & data$wave>13,])
      summary(mlm.0.l)

      mlm.1.l <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | fullint) + (1 | wave), data=data[data$pro_immig==1 & data$wave>13,])
      summary(mlm.1.l)
      anova(mlm.0.l, mlm.1.l)

      # Neutral
      mlm.0.n <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | wave), data=data[data$pro_immig==0 & data$wave>13,])
      summary(mlm.0.n)

      mlm.1.n <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | fullint) + (1 | wave), data=data[data$pro_immig==0 & data$wave>13,])
      summary(mlm.1.n)
      anova(mlm.0.n, mlm.1.n)


      # Anti-Immigration
      mlm.0.r <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | wave), data=data[data$pro_immig==-1 & data$wave>13,])
      summary(mlm.0.r)

      mlm.1.r <- blmer(agree_vig.2 ~ highnumb2 + easyacc2 + morerights2 + diffacc + diffrights + (1 | fullint) + (1 | wave), data=data[data$pro_immig==-1 & data$wave>13,])
      summary(mlm.1.r)
      anova(mlm.0.r, mlm.1.r)

#*********
# Table 2:
      xtable(rbind(anova(mlm.0, mlm.1)[,5:8],  
             anova(mlm.0.l, mlm.1.l)[,5:8], 
             anova(mlm.0.r, mlm.1.r)[,5:8]))

#*********
# Figure 3:
      par(mfrow=c(1,4), mar=c(3, 1, 1, 1))
      names <- c("Increase Numbers", "Lower Requirements", "Extend Rights", "Requirements: Nationality - Skills", "Rights: Cultural - Welfare", "", "SD: Full Policy Interactions", "SD: Waves", "SD: Residuals")

      plot(rep(1,9), 9:1, type="n", ann=F, axes=F)

      # Full Sample
      est.0 <- c(fixef(mlm.0)[-1], NA, NA,
           sqrt(summary(mlm.0)$varcor$wave[1]),
           summary(mlm.0)$sigma)

      est.1 <- c(fixef(mlm.1)[-1], NA,
           sqrt(summary(mlm.1)$varcor$fullint[1]),
           sqrt(summary(mlm.1)$varcor$wave[1]),
           summary(mlm.1)$sigma)

      plot(est.0, c(9:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 9.5), xlim=c(-.5, .5), main="Full Sample")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0)[-1]-1.96*se.fixef(mlm.0)[-1], c(9:5)-.2, fixef(mlm.0)[-1]+1.96*se.fixef(mlm.0)[-1], c(9:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      axis(2, at=9:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(9:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(9:1)+.2, pch=19, col=c(rep(1, 6), "maroon3", rep(1,2)))
      segments(fixef(mlm.1)[-1]-1.96*se.fixef(mlm.1)[-1], c(9:5)+.2, fixef(mlm.1)[-1]+1.96*se.fixef(mlm.1)[-1], c(9:5)+.2)
      text(est.1[7], 3+.2, "***",  pos=4)

      # Pro Immigration
      est.0 <- c(fixef(mlm.0.l)[-1], NA, NA,
           sqrt(summary(mlm.0.l)$varcor$wave[1]),
           summary(mlm.0.l)$sigma)

      est.1 <- c(fixef(mlm.1.l)[-1], NA,
           sqrt(summary(mlm.1.l)$varcor$fullint[1]),
           sqrt(summary(mlm.1.l)$varcor$wave[1]),
           summary(mlm.1.l)$sigma)

      plot(est.0, c(9:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 9.5), xlim=c(-.5, .5), main="Pro-Immigration")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0.l)[-1]-1.96*se.fixef(mlm.0.l)[-1], c(9:5)-.2, fixef(mlm.0.l)[-1]+1.96*se.fixef(mlm.0.l)[-1], c(9:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      #axis(2, at=9:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(9:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(9:1)+.2, pch=19, col=c(rep(1, 6), "maroon3", rep(1,2)))
      segments(fixef(mlm.1.l)[-1]-1.96*se.fixef(mlm.1.l)[-1], c(9:5)+.2, fixef(mlm.1.l)[-1]+1.96*se.fixef(mlm.1.l)[-1], c(9:5)+.2)

      # Anti-Immigration
      est.0 <- c(fixef(mlm.0.r)[-1], NA, NA,
           sqrt(summary(mlm.0.r)$varcor$wave[1]),
           summary(mlm.0.r)$sigma)

      est.1 <- c(fixef(mlm.1.r)[-1], NA,
           sqrt(summary(mlm.1.r)$varcor$fullint[1]),
           sqrt(summary(mlm.1.r)$varcor$wave[1]),
           summary(mlm.1.r)$sigma)

      plot(est.0, c(9:1)-.2, pch=19, ylab="", xlab="Est.", axes=F, ylim=c(0.5, 9.5), xlim=c(-.5, .5), main="Anti-Immigration")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(fixef(mlm.0.r)[-1]-1.96*se.fixef(mlm.0.r)[-1], c(9:5)-.2, fixef(mlm.0.r)[-1]+1.96*se.fixef(mlm.0.r)[-1], c(9:5)-.2)
      abline(v=0, lty=1, lwd=1.5, col="grey")
      #axis(2, at=9:1, names, lty=0)
      axis(1, lty=0)
      points(est.0, c(9:1)-.2, pch=19, cex=.6, col="white")
      points(est.1, c(9:1)+.2, pch=19, col=c(rep(1, 6), "maroon3", rep(1,2)))
      segments(fixef(mlm.1.r)[-1]-1.96*se.fixef(mlm.1.r)[-1], c(9:5)+.2, fixef(mlm.1.r)[-1]+1.96*se.fixef(mlm.1.r)[-1], c(9:5)+.2)
      text(est.1[7], 3+.2, "***",  pos=4)

##########################################################################
# Local Conditionality: Most and Least Preferred Policies
##########################################################################
      
      sim.mlm.1.l <- sim(mlm.1.l, 1000)
      sim.mlm.1.r <- sim(mlm.1.r, 1000)

      polpref.l <- c(as.matrix(coef(mlm.1.l)$fullint[1,1:6]) %*% c(1, 0, 0, 1, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[2,1:6]) %*% c(1, 1, 0, 1, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[3,1:6]) %*% c(1, 0, 0, 1, 1, 1),
               as.matrix(coef(mlm.1.l)$fullint[4,1:6]) %*% c(1, 1, 0, 1, 1, 1),
               
               as.matrix(coef(mlm.1.l)$fullint[5,1:6]) %*% c(1, 0, 1, 1, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[6,1:6]) %*% c(1, 1, 1, 1, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[7,1:6]) %*% c(1, 0, 1, 1, 1, 1),
               as.matrix(coef(mlm.1.l)$fullint[8,1:6]) %*% c(1, 1, 1, 1, 1, 1),
               
               as.matrix(coef(mlm.1.l)$fullint[9,1:6]) %*% c(1, 0, 0, 1, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[10,1:6]) %*% c(1, 1, 0, 1, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[11,1:6]) %*% c(1, 0, 0, 1, 1, 0),
               as.matrix(coef(mlm.1.l)$fullint[12,1:6]) %*% c(1, 1, 0, 1, 1, 0),
               
               as.matrix(coef(mlm.1.l)$fullint[13,1:6]) %*% c(1, 0, 1, 1, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[14,1:6]) %*% c(1, 1, 1, 1, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[15,1:6]) %*% c(1, 0, 1, 1, 1, 0),
               as.matrix(coef(mlm.1.l)$fullint[16,1:6]) %*% c(1, 1, 1, 1, 1, 0),
               
               as.matrix(coef(mlm.1.l)$fullint[17,1:6]) %*% c(1, 0, 0, 0, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[18,1:6]) %*% c(1, 1, 0, 0, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[19,1:6]) %*% c(1, 0, 0, 0, 1, 1),
               as.matrix(coef(mlm.1.l)$fullint[20,1:6]) %*% c(1, 1, 0, 0, 1, 1),
               
               as.matrix(coef(mlm.1.l)$fullint[21,1:6]) %*% c(1, 0, 1, 0, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[22,1:6]) %*% c(1, 1, 1, 0, 0, 1),
               as.matrix(coef(mlm.1.l)$fullint[23,1:6]) %*% c(1, 0, 1, 0, 1, 1),
               as.matrix(coef(mlm.1.l)$fullint[24,1:6]) %*% c(1, 1, 1, 0, 1, 1),
               
               as.matrix(coef(mlm.1.l)$fullint[25,1:6]) %*% c(1, 0, 0, 0, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[26,1:6]) %*% c(1, 1, 0, 0, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[27,1:6]) %*% c(1, 0, 0, 0, 1, 0),
               as.matrix(coef(mlm.1.l)$fullint[28,1:6]) %*% c(1, 1, 0, 0, 1, 0),
               
               as.matrix(coef(mlm.1.l)$fullint[29,1:6]) %*% c(1, 0, 1, 0, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[30,1:6]) %*% c(1, 1, 1, 0, 0, 0),
               as.matrix(coef(mlm.1.l)$fullint[31,1:6]) %*% c(1, 0, 1, 0, 1, 0),
               as.matrix(coef(mlm.1.l)$fullint[32,1:6]) %*% c(1, 1, 1, 0, 1, 0))

      sim.polpref.l <- cbind(coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 1, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 1,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 1, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 2,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 1, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 3,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 1, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 4,],
               
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 1, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 5,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 1, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 6,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 1, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 7,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 1, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 8,],
               
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 1, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 9,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 1, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 10,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 1, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 11,],
                   coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 1, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 12,],
               
                   coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 1, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 13,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 1, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 14,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 1, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 15,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 1, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 16,],
               
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 0, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 17,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 0, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 18,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 0, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 19,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 0, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 20,],
               
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 0, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 21,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 0, 0, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 22,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 0, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 23,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 0, 1, 1) + coef(sim.mlm.1.l)$ranef$fullint[, 24,],
               
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 0, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 25,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 0, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 26,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 0, 0, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 27,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 0, 0, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 28,],
               
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 0, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 29,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 0, 0, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 30,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 0, 1, 0, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 31,],
               coef(sim.mlm.1.l)$fixef %*% c(1, 1, 1, 0, 1, 0) + coef(sim.mlm.1.l)$ranef$fullint[, 32,])



      polpref.r <- c(as.matrix(coef(mlm.1.r)$fullint[1,1:6]) %*% c(1, 0, 0, 1, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[2,1:6]) %*% c(1, 1, 0, 1, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[3,1:6]) %*% c(1, 0, 0, 1, 1, 1),
               as.matrix(coef(mlm.1.r)$fullint[4,1:6]) %*% c(1, 1, 0, 1, 1, 1),
               
               as.matrix(coef(mlm.1.r)$fullint[5,1:6]) %*% c(1, 0, 1, 1, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[6,1:6]) %*% c(1, 1, 1, 1, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[7,1:6]) %*% c(1, 0, 1, 1, 1, 1),
               as.matrix(coef(mlm.1.r)$fullint[8,1:6]) %*% c(1, 1, 1, 1, 1, 1),
               
               as.matrix(coef(mlm.1.r)$fullint[9,1:6]) %*% c(1, 0, 0, 1, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[10,1:6]) %*% c(1, 1, 0, 1, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[11,1:6]) %*% c(1, 0, 0, 1, 1, 0),
               as.matrix(coef(mlm.1.r)$fullint[12,1:6]) %*% c(1, 1, 0, 1, 1, 0),
               
               as.matrix(coef(mlm.1.r)$fullint[13,1:6]) %*% c(1, 0, 1, 1, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[14,1:6]) %*% c(1, 1, 1, 1, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[15,1:6]) %*% c(1, 0, 1, 1, 1, 0),
               as.matrix(coef(mlm.1.r)$fullint[16,1:6]) %*% c(1, 1, 1, 1, 1, 0),
               
               as.matrix(coef(mlm.1.r)$fullint[17,1:6]) %*% c(1, 0, 0, 0, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[18,1:6]) %*% c(1, 1, 0, 0, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[19,1:6]) %*% c(1, 0, 0, 0, 1, 1),
               as.matrix(coef(mlm.1.r)$fullint[20,1:6]) %*% c(1, 1, 0, 0, 1, 1),
               
               as.matrix(coef(mlm.1.r)$fullint[21,1:6]) %*% c(1, 0, 1, 0, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[22,1:6]) %*% c(1, 1, 1, 0, 0, 1),
               as.matrix(coef(mlm.1.r)$fullint[23,1:6]) %*% c(1, 0, 1, 0, 1, 1),
               as.matrix(coef(mlm.1.r)$fullint[24,1:6]) %*% c(1, 1, 1, 0, 1, 1),
               
               as.matrix(coef(mlm.1.r)$fullint[25,1:6]) %*% c(1, 0, 0, 0, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[26,1:6]) %*% c(1, 1, 0, 0, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[27,1:6]) %*% c(1, 0, 0, 0, 1, 0),
               as.matrix(coef(mlm.1.r)$fullint[28,1:6]) %*% c(1, 1, 0, 0, 1, 0),
               
               as.matrix(coef(mlm.1.r)$fullint[29,1:6]) %*% c(1, 0, 1, 0, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[30,1:6]) %*% c(1, 1, 1, 0, 0, 0),
               as.matrix(coef(mlm.1.r)$fullint[31,1:6]) %*% c(1, 0, 1, 0, 1, 0),
               as.matrix(coef(mlm.1.r)$fullint[32,1:6]) %*% c(1, 1, 1, 0, 1, 0))


      sim.polpref.r <- cbind(coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 1, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 1,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 1, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 2,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 1, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 3,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 1, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 4,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 1, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 5,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 1, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 6,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 1, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 7,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 1, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 8,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 1, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 9,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 1, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 10,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 1, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 11,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 1, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 12,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 1, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 13,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 1, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 14,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 1, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 15,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 1, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 16,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 0, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 17,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 0, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 18,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 0, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 19,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 0, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 20,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 0, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 21,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 0, 0, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 22,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 0, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 23,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 0, 1, 1) + coef(sim.mlm.1.r)$ranef$fullint[, 24,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 0, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 25,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 0, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 26,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 0, 0, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 27,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 0, 0, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 28,],
                       
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 0, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 29,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 0, 0, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 30,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 0, 1, 0, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 31,],
                       coef(sim.mlm.1.r)$fixef %*% c(1, 1, 1, 0, 1, 0) + coef(sim.mlm.1.r)$ranef$fullint[, 32,])


      policy.labs <- c("decrease_num.increase_edu.expand_cust",   "increase_num.increase_edu.expand_cust",   "decrease_num.increase_nat.expand_cust",  
                 "increase_num.increase_nat.expand_cust",   "decrease_num.lower_edu.expand_cust",      "increase_num.lower_edu.expand_cust",    
                 "decrease_num.lower_nat.expand_cust",      "increase_num.lower_nat.expand_cust",      "decrease_num.increase_edu.expand_welf",  
                 "increase_num.increase_edu.expand_welf",   "decrease_num.increase_nat.expand_welf",   "increase_num.increase_nat.expand_welf",  
                 "decrease_num.lower_edu.expand_welf",      "increase_num.lower_edu.expand_welf",      "decrease_num.lower_nat.expand_welf",     
                 "increase_num.lower_nat.expand_welf",      "decrease_num.increase_edu.restrict_cust", "increase_num.increase_edu.restrict_cust",
                 "decrease_num.increase_nat.restrict_cust", "increase_num.increase_nat.restrict_cust", "decrease_num.lower_edu.restrict_cust",   
                 "increase_num.lower_edu.restrict_cust",    "decrease_num.lower_nat.restrict_cust",    "increase_num.lower_nat.restrict_cust",   
                 "decrease_num.increase_edu.restrict_welf", "increase_num.increase_edu.restrict_welf", "decrease_num.increase_nat.restrict_welf",
                 "increase_num.increase_nat.restrict_welf", "decrease_num.lower_edu.restrict_welf",    "increase_num.lower_edu.restrict_welf",   
                 "decrease_num.lower_nat.restrict_welf",    "increase_num.lower_nat.restrict_welf") 


      plot(polpref.r, polpref.l, xlim=c(0, 1), ylim=c(0, .5))
      text(polpref.r, polpref.l, policy.labs, pos=4, cex=.6)

      ord <- order(polpref.l)

      plot(rep(1, 32), polpref.l[ord], xlim=c(0, 3), ylim=c(0, .8))
      points(rep(2, 32), polpref.r[ord])

      segments(rep(1, 32), polpref.l[ord], rep(2, 32), polpref.r[ord])


      colnames(sim.polpref.r) <- policy.labs
      colnames(sim.polpref.l) <- policy.labs

      p.means.l <- apply(sim.polpref.l, 2, mean)
      ci.means.l <- apply(sim.polpref.l, 2, quantile, c(.025, .975))
      ord.l <- order(p.means.l, decreasing=T)

      p.means.r <- apply(sim.polpref.r, 2, mean)
      ci.means.r <- apply(sim.polpref.r, 2, quantile, c(.025, .975))
      ord.r <- order(p.means.r, decreasing=T)

#***********
# Figure A9:
      par(mar=c(3, 12, 1, 1), mfrow=c(1, 1), cex.lab=.8)
      plot(p.means.l[ord.l], c(32:1), xlim=c(0, .5), axes=F, ylab="", xlab="Policy Support", pch=19, main="Preference Order of Pro-Immigration Subset")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(ci.means.l[1, ord.l], c(32:1), ci.means.l[2, ord.l], c(32:1))
      axis(2, at=c(32:1), labels=policy.labs[ord.l], cex.axis=.8)
      axis(1)

#***********
# Figure A10:
      plot(p.means.r[ord.r], c(32:1), xlim=c(0, .8), axes=F, ylab="", xlab="Policy Support", pch=19, main="Preference Order of Anti-Immigration Subset")
      grid(lty=1, lwd=.5, col="lightgrey")
      segments(ci.means.r[1, ord.r], c(32:1), ci.means.r[2, ord.r], c(32:1))
      axis(2, at=c(32:1), labels=policy.labs[ord.r], cex.axis=.6)
      axis(1)

################################################################################################
# STUDY 3
################################################################################################

    setwd("~/YOURPATH/")
    dat <- read.dta13("replication_data_numbers_selectivity_rights_study_3.dta")
    dat <- dat[dat$attention_check_1!="",]
    dat <- dat[dat$attention_check_1=="Stimme überhaupt nicht zu",]

###############################################################################################
# Recode Variables
###############################################################################################
    
    dat$female <- ifelse(dat$gender=="Weiblich", 1, 0)
    dat$age <- 2022-dat$age
    dat$edu <- ifelse(dat$education=="Abitur, allgemeine oder fachgebundene Hochschulreife bzw. Erweiterte Oberschule der ehem. DDR mit Abschluss 12. Klasse", 1, 0)
    dat$left_right <- as.numeric(dat$left_right)

    dat$christherit[dat$cultural_lib_cons_3=="Sehr guter Vorschlag"] <- 5
    dat$christherit[dat$cultural_lib_cons_3=="Ziemlich guter Vorschlag"] <- 4
    dat$christherit[dat$cultural_lib_cons_3=="Weder guter noch schlechter Vorschlag"] <- 3
    dat$christherit[dat$cultural_lib_cons_3=="Ziemlich schlechter Vorschlag"] <- 2
    dat$christherit[dat$cultural_lib_cons_3=="Sehr schlechter Vorschlag"] <- 1

    dat$multicult[dat$cultural_lib_cons_4=="Sehr guter Vorschlag"] <- 5
    dat$multicult[dat$cultural_lib_cons_4=="Ziemlich guter Vorschlag"] <- 4
    dat$multicult[dat$cultural_lib_cons_4=="Weder guter noch schlechter Vorschlag"] <- 3
    dat$multicult[dat$cultural_lib_cons_4=="Ziemlich schlechter Vorschlag"] <- 2
    dat$multicult[dat$cultural_lib_cons_4=="Sehr schlechter Vorschlag"] <- 1


#############################################################################################
### Tradeoff A
#############################################################################################
    
    dat$tradeoff_a <- ifelse(dat$tradeoff_a=="" | dat$tradeoff_a=="NA" | dat$tradeoff_a=="Weiß nicht", NA, dat$tradeoff_a)
    dat$tradeoff_a[dat$tradeoff_a=="Voll und ganz unterstützen"] <- 6
    dat$tradeoff_a[dat$tradeoff_a=="sehr unterstützen"] <- 5
    dat$tradeoff_a[dat$tradeoff_a=="etwas unterstützen"] <- 4
    dat$tradeoff_a[dat$tradeoff_a=="etwas ablehnen"] <- 3
    dat$tradeoff_a[dat$tradeoff_a=="sehr ablehnen"] <- 2
    dat$tradeoff_a[dat$tradeoff_a=="voll und ganz ablehnen"] <- 1
    dat$tradeoff_a <- as.numeric(dat$tradeoff_a)

    dat$tradeoff_a1 <- ifelse(dat$tradeoff_a1=="dem Mittleren Osten und Nordafrika", 1, 0)
    dat$tradeoff_a2 <- ifelse(dat$tradeoff_a2=="erhöhen", 1, 0)

#############################################################################################
# Balance Check
#############################################################################################
    
#**********    
# Table A7:
      summary(lm(female ~ tradeoff_a1*tradeoff_a2, data=dat))
      summary(lm(age ~ tradeoff_a1*tradeoff_a2, data=dat))
      summary(lm(edu ~ tradeoff_a1*tradeoff_a2, data=dat))
      summary(lm(left_right ~ tradeoff_a1*tradeoff_a2, data=dat))
      summary(lm(christherit ~ tradeoff_a1*tradeoff_a2, data=dat))
      summary(lm(multicult ~ tradeoff_a1*tradeoff_a2, data=dat))

      table(interaction(dat$tradeoff_a1, dat$tradeoff_a2))


#********** 
# Table A4:
      stargazer(dat[,c("female", "age", "edu", "left_right", "christherit", "multicult")])


      summary(lm(tradeoff_a ~ tradeoff_a1 + tradeoff_a2, data=dat))
      summary(lm(tradeoff_a ~ tradeoff_a1*tradeoff_a2, data=dat))
      aggregate(dat$tradeoff_a, by=list(dat$tradeoff_a1, dat$tradeoff_a2), mean, na.rm=T)


      dat$tradeoff_a.2 <- ifelse(dat$tradeoff_a<=3, 0, 1)

      prop.table(table(dat$tradeoff_a2, dat$tradeoff_a.2))
      round(prop.table(table(dat$tradeoff_a2, dat$tradeoff_a.2)), 3)

      dat$genpref[dat$tradeoff_a2==0 & dat$tradeoff_a.2==0] <- "oppose_decrease"
      dat$genpref[dat$tradeoff_a2==1 & dat$tradeoff_a.2==0] <- "oppose_increase"
      dat$genpref[dat$tradeoff_a2==0 & dat$tradeoff_a.2==1] <- "support_decrease"
      dat$genpref[dat$tradeoff_a2==1 & dat$tradeoff_a.2==1] <- "support_increase"

      
##############################################################################################
### Tradeoff B
##############################################################################################
      
      dat$tradeoff_b <- ifelse(dat$tradeoff_b=="" | dat$tradeoff_b=="NA" | dat$tradeoff_b=="Weiß nicht", NA, dat$tradeoff_b)
      dat$tradeoff_b[dat$tradeoff_b=="Voll und ganz unterstützen"] <- 6
      dat$tradeoff_b[dat$tradeoff_b=="sehr unterstützen"] <- 5
      dat$tradeoff_b[dat$tradeoff_b=="etwas unterstützen"] <- 4
      dat$tradeoff_b[dat$tradeoff_b=="etwas ablehnen"] <- 3
      dat$tradeoff_b[dat$tradeoff_b=="sehr ablehnen"] <- 2
      dat$tradeoff_b[dat$tradeoff_b=="voll und ganz ablehnen"] <- 1
      dat$tradeoff_b <- as.numeric(dat$tradeoff_b)

      dat$tradeoff_changepref <- dat$tradeoff_b - dat$tradeoff_a

      round(mean(dat$tradeoff_changepref>0, na.rm=T), 2)
      round(mean(dat$tradeoff_changepref<0, na.rm=T), 2)
      round(mean(dat$tradeoff_changepref==0, na.rm=T), 2)

      sum(dat$tradeoff_changepref>0, na.rm=T)
      sum(dat$tradeoff_changepref<0, na.rm=T)
      sum(dat$tradeoff_changepref==0, na.rm=T)

      dat$change[dat$tradeoff_changepref>0] <- 1
      dat$change[dat$tradeoff_changepref<0] <- -1
      dat$change[dat$tradeoff_changepref==0] <- 0
      table(dat$change)
      prop.table(table(dat$change))


      table(dat$genpref, dat$change)


      dat$tradeoff_b1 <- ifelse(dat$tradeoff_b1=="", NA, dat$tradeoff_b1)
      dat$tradeoff_b2 <- ifelse(dat$tradeoff_b2=="", NA, dat$tradeoff_b2)
      dat$tradeoff_b3 <- ifelse(dat$tradeoff_b3=="", NA, dat$tradeoff_b3)
      dat$tradeoff_b4 <- ifelse(dat$tradeoff_b4=="", NA, dat$tradeoff_b4)
      dat$tradeoff_b5 <- ifelse(dat$tradeoff_b5=="", NA, dat$tradeoff_b5)
      dat$tradeoff_b6 <- ifelse(dat$tradeoff_b6=="", NA, dat$tradeoff_b6)
      dat$tradeoff_b7 <- ifelse(dat$tradeoff_b7=="", NA, dat$tradeoff_b7)

###########################################################################################
# Balance Check
###########################################################################################
      
#**********      
# Table A8:
        summary(lm(female ~ tradeoff_b5*tradeoff_b6, data=dat))
        summary(lm(age ~ tradeoff_b5*tradeoff_b6, data=dat))
        summary(lm(edu ~ tradeoff_b5*tradeoff_b6, data=dat))
        summary(lm(left_right ~ tradeoff_b5*tradeoff_b6, data=dat))
        summary(lm(christherit ~ tradeoff_b5*tradeoff_b6, data=dat))
        summary(lm(multicult ~ tradeoff_b5*tradeoff_b6, data=dat))


        dat$more_selective <- ifelse(dat$tradeoff_b5=="verschärfen", 1, 0)
        dat$less_selective <- ifelse(dat$tradeoff_b5=="erleichtern", 1, 0)

        dat$less_rights <- ifelse(dat$tradeoff_b6=="einschränken", 1, 0)
        dat$more_rights <- ifelse(dat$tradeoff_b6=="ausweiten", 1, 0)

        dat$pref_increase <- ifelse(dat$tradeoff_b1=="einverstanden" & dat$tradeoff_b2=="erhöht"
                           | dat$tradeoff_b1=="nicht einverstanden" & dat$tradeoff_b2=="verringert", 1, 0)

        dat$pref_decrease <- ifelse(dat$tradeoff_b1=="einverstanden" & dat$tradeoff_b2=="verringert"
                            | dat$tradeoff_b1=="nicht einverstanden" & dat$tradeoff_b2=="erhöht", 1, 0)

        dat$no_increase <- ifelse(dat$tradeoff_b1=="nicht einverstanden" & dat$tradeoff_b2=="erhöht", 1, 0)
        dat$no_decrease <- ifelse(dat$tradeoff_b1=="nicht einverstanden" & dat$tradeoff_b2=="verringert", 1, 0)

        chisq.test(table(dat$pref_increase, dat$change))

###########################################################################################
# Preferences
###########################################################################################
        
    # All Immigrants
    m.3.a <- (lm(tradeoff_b ~ more_selective + less_rights, data=dat[dat$pref_decrease==1,]))
    m.3.a.i <- (lm(tradeoff_b ~ more_selective*less_rights, data=dat[dat$pref_decrease==1,]))

    m.3.b <- (lm(tradeoff_changepref ~ more_selective + less_rights, data=dat[dat$pref_decrease==1,]))
    m.3.b.i <- (lm(tradeoff_changepref ~ more_selective*less_rights, data=dat[dat$pref_decrease==1,]))

#**********
# Table A11:
      stargazer(m.3.a, m.3.a.i, m.3.b, m.3.b.i)

      m.4.a <- (lm(tradeoff_b ~ less_selective + more_rights, data=dat[dat$pref_increase==1,]))
      m.4.a.i <- (lm(tradeoff_b ~ less_selective*more_rights, data=dat[dat$pref_increase==1,]))

      m.4.b <- (lm(tradeoff_changepref ~ less_selective + more_rights, data=dat[dat$pref_increase==1,]))
      m.4.b.i <- (lm(tradeoff_changepref ~ less_selective*more_rights, data=dat[dat$pref_increase==1,]))

#**********
# Table A12:
      stargazer(m.4.a, m.4.a.i, m.4.b, m.4.b.i)


      # MENA
      summary(lm(tradeoff_b ~ more_selective + less_rights, data=dat[dat$pref_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      summary(lm(tradeoff_changepref ~ more_selective + less_rights, data=dat[dat$pref_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

      summary(lm(tradeoff_b ~ less_selective + more_rights, data=dat[dat$pref_increase==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      summary(lm(tradeoff_changepref ~ less_selective + more_rights, data=dat[dat$pref_increas==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

###########################################################################################
# Those Denied Their Preferences
###########################################################################################
      
      # All Immigrants
      m.1.a <- (lm(tradeoff_b ~ more_selective + less_rights, data=dat[dat$no_increase==1,]))
      m.1.a.i <- (lm(tradeoff_b ~ more_selective*less_rights, data=dat[dat$no_increase==1,]))

      m.1.b <- (lm(tradeoff_changepref ~ more_selective + less_rights, data=dat[dat$no_increase==1,]))
      m.1.b.i <- (lm(tradeoff_changepref ~ more_selective*less_rights, data=dat[dat$no_increase==1,]))

#***********
# Table 3:
      stargazer(m.1.a, m.1.a.i, m.1.b, m.1.b.i)


      # MENA
      m.1.c <- (lm(tradeoff_b ~ more_selective + less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      m.1.c.i <- (lm(tradeoff_b ~ more_selective*less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

      m.1.d <- (lm(tradeoff_changepref ~ more_selective + less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      m.1.d.i <- (lm(tradeoff_changepref ~ more_selective*less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

#***********
# Table A13:
      stargazer(m.1.c, m.1.c.i, m.1.d, m.1.d.i)

      # Europe
      m.1.e <- (lm(tradeoff_b ~ more_selective + less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))
      m.1.e.i <- (lm(tradeoff_b ~ more_selective*less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))

      m.1.f <- (lm(tradeoff_changepref ~ more_selective + less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))
      m.1.f.i <- (lm(tradeoff_changepref ~ more_selective*less_rights, data=dat[dat$no_increase==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))

#***********
# Table A14:
      stargazer(m.1.e, m.1.e.i, m.1.f, m.1.f.i)


      # All Immigrants
      m.2.a <- (lm(tradeoff_b ~ less_selective + more_rights, data=dat[dat$no_decrease==1,]))
      m.2.a.i <- (lm(tradeoff_b ~ less_selective*more_rights, data=dat[dat$no_decrease==1,]))

      m.2.b <- (lm(tradeoff_changepref ~ less_selective + more_rights, data=dat[dat$no_decrease==1,]))
      m.2.b.i <- (lm(tradeoff_changepref ~ less_selective*more_rights, data=dat[dat$no_decrease==1,]))

#***********
# Table 4:
      stargazer(m.2.a, m.2.a.i, m.2.b, m.2.b.i)

      # MENA
      m.2.c <- (lm(tradeoff_b ~ less_selective + more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      m.2.c.i <- (lm(tradeoff_b ~ less_selective*more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

      m.2.d <- (lm(tradeoff_changepref ~ less_selective + more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))
      m.2.d.i <- (lm(tradeoff_changepref ~ less_selective*more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3=="dem Mittleren Osten und Nordafrika",]))

#***********
# Table A15:
      stargazer(m.2.c, m.2.c.i, m.2.d, m.2.d.i)

      # Europe
      m.2.e <- (lm(tradeoff_b ~ less_selective + more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))
      m.2.e.i <- (lm(tradeoff_b ~ less_selective*more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))

      m.2.f <- (lm(tradeoff_changepref ~ less_selective + more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))
      m.2.f.i <- (lm(tradeoff_changepref ~ less_selective*more_rights, data=dat[dat$no_decrease==1 & dat$tradeoff_b3!="dem Mittleren Osten und Nordafrika",]))

#***********
# Table A16:
      stargazer(m.2.e, m.2.e.i, m.2.f, m.2.f.i)

################################################################################################



